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ABSTRACT 

QSO near-zones are an important probe of the the ionization state of the IGM at 
2 ; ^ 6 — 7, at the end of reionization. We present here high-resolution cosmological 
3D radiative transfer simulations of QSO environments for a wide range of host halo 
masses, Our simulated near-zones reproduce both the overall decrease 

of observed near-zone sizes at 6 < z < 7 and their scatter. The observable near-zone 
properties in our simulations depend only very weakly on the mass of the host halo. The 
size of the H ll region expanding into the IGM is generally limited by (super-)Lyman 
Limit systems loosely associated with (low-mass) dark matter haloes. This leads to a 
strong dependence of near-zone size on direction and drives the large observed scatter. 
In the simulation centred on our most massive host halo, many sightlines show strong 
red damping wings even for initial volume averaged neutral hydrogen fractions as low 
as ^ 10“^. For QSO lifetimes long enough to allow growth of the central supermassive 
black hole while optically bright, we can reproduce the observed near-zone of ULAS 
J1120+0641 only with an IGM that is initially neutral. Our results suggest that larger 
samples of z > 7 QSOs will provide important constraints on the evolution of the 
neutral hydrogen fraction and thus on how late reionization ends. 

Key words: galaxies: high-redshift - quasars: absorption lines - intergalactic medium 
- dark ages, reionization, hrst stars 


1 INTRODUCTION 


An important unsolved problem in modern cosmology is un¬ 
derstanding how our Universe transitioned from the “dark 
ages”, following recombination, to the ionized Universe we 
observe today. Despite much progress in the last decade 
there is still considerable uncertainty with regard to the ex¬ 
act timing and topology of reionization and which sources 
are the main contributors to the ionizing emissivity. The evo¬ 
lution of the Lyman-a (Lya) forest and detection of a Gunn- 
Peterson trough blueward of Lya suggest that the neutral 
fraction of the intergalactic medium (IGM) is rapidly in¬ 
creasing for > 5.7 ( Fan et al.|2006 Becker et al.|2015 ). Ap¬ 
parent drops in the number of Lyman-a emitters at z = 6—7 
also suggest a rapidly increasing neutral fraction of hydro¬ 
gen (e.g., Stark et al. 2010| Ouchi et al. 2010| Finkelstein 
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et al.|[MT3 Pentericci et al.|[2014 ). Interestingly, the latest 
measurements of the optical depth due to Thomson scatter¬ 
ing from observations of the cosmic microwave background 
(CMB) are signihcantly lower than earlier measurements 
( Planck Collaboration et al.||2014| 2015). For the rather ar¬ 
tificial assumption of instantaneous reionization, the new 
inferred reionization redshift from the latest CMB measure¬ 
ments is j^reion = For morc realistic extended evolu¬ 

tionary histories of the neutral hydrogen fraction, the new 
measurements suggest that reionization ends somewhat later 
than previously thought ( Haardt fc Madau|2012 [Choudhury 


et al. 12014) [Robertson et al.|2015|(Ba 

et al.|2015|). 


et al.|2015||Chardm 


With a number of near-infrared surveys under way to 
further push the redshift barrier beyond z > 7, this makes 
QSO near-zones a promising tool to further investigate how 
reionization ended. Near-zones are defined as regions of 
transmitted flux blueward of the Lya emission, extending 
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from the source to the point where the flux drops below 10 
per cent for the first time ( Fan et al.|2006 ). UKIDSS, VISTA 
and Pan-STARRS are now returning the first discoveries of 
2 : > 6.5 quasars (Mortlock et al. [Venemans et al.|2013| 


2015). The near-zones of these bright QSOs are a unique 
observational tool to study the high-redshift IGM, as they 
allow us to push past the redshift at which the opacity in the 
Lya forest becomes too large to extract detailed information 
(e.g .,|Fan et al.||2006| |Carilli et al. poTol i. 

Carilli et al.| ( |20io| l found evidence for decreasing near¬ 


zone size with increasing redshift in the range 5.7 < z < 6.4, 
suggesting evolution in the neutral fraction of the IGM, as 
well as substantial scatter in the observed near-zone sizes. 
This trend continues to the current highest redshift QSO, 
ULAS J1120+0641, at z = 7.085 (see Venemans et al. 


(2015); Wu et al. (2015) and Reed et al. (2015) for more 


recent updates with additional QSOs at z > 6). 

The small near-zone size of ULAS J1120+0641, together 
with the possible discovery of a damping wing observed red- 
ward of Lya transmission, also suggests that the IGM at 
z ~ 7 may be significantly more neutral than at z ~ 6 
([Mortlock et al.||2011 Bolton et al.||20lT ), but see [Bosman 


& Becker (2015) for a reassessment of the significance of the 


evidence for suppressed Lyman-a emission in this QSO by 
a significantly neutral surrounding IGM. 

The relation of the sizes of near-zones to the neutral 
fraction of the surrounding IGM is, however, not straight for¬ 
ward. The size of the near-zone generally underestimates the 


size of the H ii region (e.g., Bolton & Haehnelt 2007a Maselli 
et al.|2007 ). Maselli, Ferrara V Gallerani (2009) showed that 


in a sample of z ~ 6 quasars, all the observed near-zones 
sizes appear to be due to the proximity zone around the 
quasar, rather than tracing the H ll region. Reionization is 
of course not expected to be a spatially homogeneous pro¬ 


cess ( Miralda-Escude, Haehnelt Rees]|2000 Furlanetto & 


Oh||2005| ) and, as quasars may form in biased, high-density 
regions, the surrounding IGM may already be substantially 
ionized by the time the quasar turns on ( Lidz et al.|2007 ). 

It is also not clear to what extent the environment 
and, in particular, the mass of the host halo influence the 


observed properties of high-redshift QSO near-zones ( Lidz 
et al. 2007| Wyithe, Bolton V Haehnelt 2008). Unfortunately 


observational constraints on the environment and host halo 
masses of high-redshift QSOs so far have remained incon¬ 


clusive. Bahados et al. (2013) and Simpson et al. (2014) 


find no evidence for the overdensities that would be ex¬ 
pected around the most massive haloes, while [Morselli et al] 
(2014) claim to observe an overdense environment around 
four z ^ 6 QSOs. The mass of the black hole powering 
ULAS J1120+0641 is estimated to be (2.0±J;7) x 10^ M© 
( Mortlock et al.]|2011 ). Based on high-resolution numerical 
simulations of the growth of galaxies and black holes in a 
cosmological context, it has been shown that it is possible 
to form a black hole of this mass by growth limited to the 
Eddington accretion rate from a seed of 10^ M© at z ~ 15 
but these black holes will only be found in the most mas¬ 
sive dark matter haloes (with Mhaio ~ 10^^ Mq) at z ~ 6 


( Sijacki, Springel V Haehnelt|2Q09 Gosta et al.|2014 ). Gon- 


versely, Eanidakis et al.j ( |2013| ) use semi-analytic modelling 
allowing for super-Eddington accretion rates to suggest that 
luminous quasars are found in less massive haloes; however, 
their modelling fails to reproduce the bright end of the high- 


redshift quasar luminosity function, as 10^ M© black holes 
are not formed until z < 4. We therefore simulate here three 
different regions of the IGM at z = 7, centred on host haloes 
spanning a wide range of masses from 10^° to 10^^ M© and 
covering a range of environments from an average density 
region to a very highly biased proto-cluster region. 

This work follows on from the work of IBolton et al.l 


(2011) but differs in two important respects: first, we test 


a range of host halo masses and environments extending 
to a rather large host halo mass and highly biased envi¬ 
ronments predicted by realistic models for the growth of 
billion solar mass black holes believed to power bright high- 
redshift QSOs; second, we use the 3D radiative transfer code 
RADAMESH ( Gantalupo Porciani|2()TT ) which allows us to 
sample the near-zones effectively in all directions. By choos¬ 
ing different background photo-ionization rates, we vary the 
initial neutral fraction of the gas in the surrounding IGM 
before the QSO turns on. 

This paper is organised as follows. In section 2, we de¬ 
scribe the numerical simulations we used to model the ion¬ 
ization fronts around high-redshift QSOs. In section 3, we in¬ 
vestigate the effect of different environments and initial neu¬ 
tral fractions on the QSO near-zones. In section 4, we com¬ 
pare our simulations to the observed near-zones at 2 ; = 6 — 7 
and consider how best to constrain the neutral fraction at 
2 : > 6.5. Einally, in section 5, we present our conclusions. 
We use the cosmological parameters h = 0.73, Um = 0.25, 
LIa = 0.75 and Llbh^ = 0.024. The prefixes ‘p’ and ‘c’ before 
distance measures refer to proper and comoving distances 
respectively. 


2 SIMULATIONS OF QSO NEAR-ZONES 

2.1 Hydrodynamical Simulations of High-Redshift 
Haloes 

As the mass of the dark matter haloes that bright QSOs 
lie in at 2 ; ~ 6 — 7 is still uncertain, we have simulated a 
range of cosmological environments, following the approach 

centred on a 


m 

10 


Gosta et al.j (|2014|): an 

7 


average region. 


h~ Mq halo, an “intermediate” region, centred on a 
1.7 X 10^^h~^ Mq halo and an “overdense” region, centred 
on a 2.5 x 10^^h~^ Mq halo. To model these three distinct 
regions, we resimulated three haloes taken from a snapshot 
of the Millennium simulation ( Springel et al.|20Q5 ) at 2 : ~ 6, 
using the TreePM-SPH code p-GADGEt3 (last described in 
Springel|2005 ). The Millennium simulation is a dark-matter 


only simulation of a cube with side length 500/i“^ Mpc run 
from redshift 127 up to the present time. The large volume 
of the simulation allows for formation of rare objects, such 
as the most massive halo we simulate. We used the group 
catalogues to select haloes of desired mass at 2 ; = 6 and 
traced all particles in a radius around these haloes back to 
2 ; = 127 to generate our initial conditions. We included a fac¬ 
tor of 4^ times more particles than the parent simulation in 
a high-resolution region, which was also populated with gas 
particles. Outside this region, increasingly lower resolution 
dark matter particles were also simulated so as to correctly 
capture the large-scale tidal forces while minimising com¬ 
putational cost. The spatial resolution of our simulation is 
a factor 3 lower than previous studies of 2 ; ~ 7 near-zones 
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Figure 1. Maps of the projected overdensity for the three simulated regions at 2 ; = 7: average (left), intermediate (middle) and overdense 
(right). The average region contains some thin filaments as well as large voids, while the intermediate region contains more pronounced 
filaments and also smaller voids. The overdense region contains several prominent filaments surrounding the central massive halo. The 
top panel shows the 10 h~^ cMpc surrounding the host halo, the middle panel shows a slice through the entire region and the bottom 
panel shows the projected overdensity with the dark matter haloes contained in a slice of the same thickness overlaid. 


by Bolton et al. (2011). Unfortunately, this is an unavoid¬ 
able compromise that had to made in order to run large 
boxes containing the rare, massive haloes in which we are 
interested. 

We model star formation in the same way as 


Haehnelt & Springel (2004), such that gas particles with 


density greater than 1000 times the mean baryon density 


and a temperature less than 10^ K are turned into collision¬ 
less star particles. Simulations with this prescription have 
been shown to reproduce well the low-density IGM respon¬ 
sible for the bulk of the absorption in the high-z Lyman- 
a forest, while providing increased numerical efficiency. We 
have also investigated a simulation with a more physically 
motivated star formation prescription that also included the 
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Figure 2. Left: Overdensity within a sphere with radius specified on the horizontal axis, centred on the host halo, for the three regions. 
Outside the central 2 pMpc, the intermediate and overdense regions show similar densities that are still slightly above the average cosmic 
density, while the density of the average region lies much closer to the cosmic density ai z = 7. Beyond a radius of 5 pMpc, all three 
regions display similar overdensities. Middle: The number of haloes above a cut-off mass log MHalo(^~^^o) = 8-^ (solid line) and 10 
(dotted line) contained in a sphere with radius specified on the horizontal axis for the three simulations. Right: The median neutral 
fraction over all sight lines for our three regions with initial volume-weighted average neutral fraction {/h i}init =0.02 after 1 Myr (top 
panel) and 10 Myr (bottom panel). 


Region 

MHalo(xlOl2fe-lM0) 

Ar=0.5 

Ar=3 

Average 

0.01 

1.00 

0.95 

Intermediate 

0.17 

1.63 

1.11 

Overdense 

2.50 

2.86 

1.10 


Table 1. Here we list the mass of the host halo and the overden¬ 
sity A = p/p (where p is the critical density) of a sphere with 
radius 0.5 and 3 pMpc, centred on the host halo, for the three 
regions we simulate at 2 ; = 7. 


effect of stellar and AGN feedback, which we discuss in Ap¬ 
pendix 

The radius of our high-resolution region is approxi¬ 
mately 40 h~^ comoving Mpc (the exact radius of each simu¬ 
lation depends on the virial radius of the halo of interest) and 
the mass of the gas particles is 2.8 x 10® h~^ M©. The soft¬ 
ening length is 1.25 h~^ comoving kpc. The haloes contain 
at least 32 dark matter particles and the virial temperature 
of the smallest haloes is Ty ~ 650 K for /i = 0.61. 

The Haardt & Madau (2012) model of the UV back¬ 
ground was used, assuming the optically thin limit for ion¬ 
izing radiation. Note that although we will later recompute 
the ionization state of the gas for a series of different photo¬ 
ionization rates, the change of UV background should in 
principle affect the temperature of the gas, an effect which 
we neglect when we reprocess the simulations. 

Examples of the density field in the immediate vicin¬ 
ity of the three host haloes and in a slice through the full 
simulation volume are shown in the top and middle panels 
of Figure The average density region is characterised by 
thinner filaments of gas and large voids, while the overdense 
region shows more prominent filaments and higher density 
diffuse gas. Within the central 2 pMpc, the difference in the 
mean overdensity of the three regions is evident (see also 
the left panel of Figure]^, with the mean overdensity of the 
region increasing with increasing mass of the host halo as 
expected. Beyond this distance, the intermediate and over- 
dense regions display similar mean densities, while the mean 


density of the average region remains lower until a radius of 
5 pMpc. This is described more quantitatively in Table ^ 
In the bottom panel of Figurewe plot the distribution 
of the subhaloes in a slice with thickness 0.03 pMpc. Note 
that the most massive halo in the average region is not our 
host halo (see, for example, the 8 x 10^® h~^ M© halo in the 
lower right of the bottom left image in Figure [^. However, 
in the region surrounding the host halo, the distribution of 
haloes is clearly different (middle panel of Figure]^. We find 
that the regions centred on more massive haloes also contain 
more haloes of smaller masses. This is true for both haloes 
with masses above log Mq) = 8.5 and 10. Similar 

to the trend with overdensity in the left panel of Figure 
we find that the number of haloes in the region is increasing 
with the mass of the host halo up to a radius of about 2 
pMpc, after which the intermediate and overdense regions 
become very alike. 


2.2 Modelling QSO Ionization Fronts 


To calculate the neutral fraction of the gas after it is ionized 
by the QSO, we post-processed the simulations with the 
3D radiative transfer code RADAMESH. RADAMESH (Radia¬ 
tive transfer on an ADAptive MESH) is a three-dimensional 
radiative transfer code, described in detail in |Cantalupo| 
& Porciani (2011). It makes use of a photon-conserving. 


ray-tracing approach that is adaptive in space and time. 
This grid-based algorithm combines elements from two com¬ 
mon techniques for solving radiative transfer problems: the 
Monte Carlo approach and the method of characteristics. 
This scheme allows for accurate and computationally effi¬ 
cient solutions of the radiative transfer equation. By track¬ 
ing the path of ionizing radiation over a static density field, 
the time evolution of the neutral hydrogen and helium frac¬ 
tions, as well as the temperature, can be calculated. It has 
been tested on the standard tests for radiative transfer codes 


from Iliev et al. (2006) and shown to perform well. 


Fully coupling simulations of radiative transfer to the 
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hydro dynamical simulations is a computationally expensive 
problem, but fortunately for the problem considered here, 
this is of little relevance. For the expansion of a H ll re¬ 
gion around a luminous QSO, the relevant time scale for 
the motion of the ionization front is much shorter than the 
hydrodynamical time scales for the gas. We have also inves¬ 
tigated the timescale for photoevaporation of gas in haloes 
with masses equivalent to the minimum mass our simulation 
resolves using the method outlined in |Iliev, Shapiro fc Raga| 
(2005), and find that the time for evaporation is on the order 
of 100 Myr for a halo with Mnaio = lO^/i^^M© at a radius 
10 comoving Mpc from the source. As the time interval over 
which we evolve the radiative transfer is much shorter than 
this, we do not expect photoevaporation to be an issue. 

For simplicity and for computational speed, we model 
the initial ionization state of the IGM using a homoge¬ 
neous UV background together with a self-shielding model 
to account for optically thick absorbers. The effect of self¬ 
shielding on the initial ionization state was implemented in 
post-processing using a simple self-shielding density thresh¬ 
old implemented similarly to Rahmati et al. (2013). Al¬ 
though fluctuations in the UV background are expected to 
be smaller at low redshift, when the mean free path of an 


ionizing photon is short, as expected for z > 6 (e.g., Songaila 


& Cowie 2010), these fluctuations can become large and 


modelling observations with a uniform UV background can 
become difficult (Becker et al. 2015 Chardin et al. 2015). 
It would obviously be more realistic when considering the 
ionization state of the IGM to take into account the con¬ 
tribution to the ionizing background of individual sources, 
however to model this correctly would add an additional 
level of complexity to this problem and is beyond the scope 
of this paper. 

To perform the radiative transfer, outputs from the hy¬ 
drodynamical simulations at z = 6 and 7 were divided into 
octants and mapped onto multi-mesh grids, each with base 
grid size 512^. The grid was refined wherever 16 or more 
gas particles were present in a cell, resulting in four levels 
of refinement. Cells on the highest level of refinement had 
a spatial resolution of 0.95 pkpc. The initial temperature 
of the gas was taken from the hydrodynamical simulation. 
Collisional ionizations are also taken into account. Taking 
the same values for the properties of ULAS J1120+0641 as 


Bolton et al. (2011), we placed a source with ionizing photon 


densest cell of the central halo. By keeping the properties of 
the QSO fixed and varying the mass of the host halo and the 
redshift, we were able to estimate the effect that the mass of 
the host halo and the evolution of the gas distribution has 
on the extent of the near-zones. 

We consider a range of initial ionization states for the 
gas, giving volume-weighted initial neutral fractions rang¬ 
ing from (/h i)init ~ 10“^ to 1 for the neutral fraction in 
the vicinity of the quasar. The photo-ionization rates and 
corresponding volume-weighted neutral fractions are sum¬ 
marised in Table 12] Photo-ionization rates for He I and He ll 
were calculated by rescaling to match the ratios between the 
photo-ionization rates in Haardt V Madau (2012). 

Ionizing UV radiation was sampled from twenty loga¬ 
rithmically spaced energy bins, ranging from 1 to 8 Ryd. A 
time-dependent, non-equilibrium chemistry solver computes 
the number density of six species (H I, H ll. He I, He ll. 


logrHi(s Q ifuiiz = 6))init ifuiiz = 7))init 


-17.0 

9.78 X 10-1 

9.82 X 10-1 

-14.8 

5.21 X 10-2 

1.11 X 10-1 

-14.3 

9.42 X 10-3 

1.99 X 10-2 

-13.7 

1.43 X 10-3 

2.20 X 10-3 

-12.8 

1.14 X 10-"! 

1.68 X 10-4 


Table 2. H i photo-ionization rates used to calculate the initial 
ionization state of the gas with the [Rahmati et al.| ( |2013| ) model 
for self-shielded absorbers at 2 ; = 6 and 7. The volume weighted 
neutral hydrogen fraction as calculated in the overdense region is 
also shown. 


He III, e~) as well as the temperature of the gas. Recombi¬ 
nation radiation was treated using an “on-the-spot” approx¬ 
imation, whereby ionizing photons are absorbed in the same 
cell in which they were emitted (i.e., that the mean free path 
of the ionizing photon is smaller than the minimum length 
scale resolved by our simulation) and that the photons are 
absorbed by the same species from which they were emitted. 

We further assume the speed of light in our simulation 
to be infinite. Although this means that the speed that the 
H II region initially expands at can appear to exceed the 
speed of light, it can be shown that it is a good description 
of how the front propagates from the point of view of an 
observer for spectra along the li ne of sight ([White et ah 
2003] [Yu[ 2005[ [Shapiro et al.[2006[ [Bolton V Haehnelt[2007a[ 


Gantalupo, Porciani V Lilly 2008). The QSO lifetime tq 


that we will refer to throughout this paper is the time at 
which the photons were emitted from the source. This can 
be related to the time t' by tq = — R{t')/c, where R{t') 

is the distance that has been travelled in a time t' and c is 
the speed of light. 


3 QSO NEAR-ZONE SIZES 
3.1 Evolution of the H II Regions 


In the top three panels of Figure]^ we show the neutral hy¬ 
drogen fraction of the gas in the three regions for three of 
the initial average neutral fractions we consider after tq = 1 
Myr at z = 7. In the cases with (/h i)init = 0.002 and 0.02 
(top and second panels), there is a large scatter along differ¬ 
ent sightlines in the distance the ionization front has trav¬ 
elled, with expansion preferentially into the less dense re¬ 
gions. This is in contrast to the ionization front modelled by 


Maselli et al. (2007) at z = 6.1, whose simulated H ii region 


does not show strong deviations from a spherical ionization 
front. This is perhaps because our higher-resolution simula¬ 
tions are resolving smaller, dense clumps of neutral gas that 
produce a shadowing effect, although this effect may become 
less pronounced if we included the effect of scattered photons 
from recombinations. In the third panel of Figurewe show 
the case with (/h i)init = 1.0. Here we find a more isotropic 
expansion of the ionization front into the surrounding IGM, 
with little scatter in the radius of the ionization front and a 
generally smaller H ii region. After 5 Myr (bottom panel), 
although the size of the H ii region and its scatter have in¬ 
creased, the ionized region is still small compared to the top 
two panels. 

The three columns of Figure]^ show our average (left), 
intermediate (middle) and overdense (right) regions. For 
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Figure 3. Maps of the projected neutral fraction for average (left), intermediate (middle) and overdense (right) regions over a slice 0.03 
pMpc thick at 2 ; = 7. The top row shows simulations after Iq = 1 Myr with initial volume-weighted neutral fraction (/h i}init = 0.002, 
second row has {/h i}init = 0-02 and third row has {/h i}init = l-O- The fourth row has {/h i}init = 1-0 and is shown at tg = 5 Myr. Note 
that this is not an instantaneous picture, but rather shows the size of the H ll region that would be measured by an observer along each 
sightline at fixed quasar lifetime. 
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Figure 4. Left: The projected neutral fraction over a slice with thickness 0.03 pMpc for the overdense region with initial neutral fraction 
(/h i}init = 0.02. (Sub)haloes in a slice of the same thickness found using subfind are overplotted. Middle: H l column density at the 
position where the near-zone is found over a range ±500 km s“^ after Iq = 1 Myr in our overdense region with {/h i}init = 0.02. The 
different colours represent the average neutral fraction in the interval we integrate over and the values are represented in the colour 
bar in the right panel. Right: As before, but for total hydrogen column density. The dashed black line indicates the column density 
corresponding to the mean density that we would expect in this velocity interval. 





log ' Mq) log ' M^) log N„ (cm'^) 


Figure 5. Left: A histogram of the masses of the halo closest to the edge of the ionization front after Iq = 1 Myr in our overdense region 
(black line) and the masses of the dark matter haloes found in our simulation (red line). We find that the lower mass haloes tend to be 
found nearest to the ionization front. Middle: The distance from the edge of the ionization front to the nearest halo, d, normalised by the 
virial radius, R 20 O 5 against the mass of that halo. The blue line is the median value. The red line is the median value of data measuring 
the distance to the nearest halo for a random point in our simulation. The contribution from the host halo is not seen because, in that 
case, d = 0. Right: The distance from the edge of the front to the nearest halo d, normalised by the virial radius, R 20 O 5 against the H l 
column density in a ±500 km s“^ range around the edge of the ionization front after Iq = 1 Myr in our overdense region. The points 
are coloured by the total hydrogen column density. 


each of the initial ionization states, the three regions display 
a similar maximum distance that the front has travelled and 
also a similar scatter. This is also shown in the right panel 
of Figure where we plot the median neutral fraction of 
all sightlines for our three density fields after the quasar has 
been on for 1 Myr in the simulation with initial volume- 
weighted neutral fraction (/Hi)init = 0.02. Up to a neutral 
fraction log /hi ~ —5, there is little difference between the 
three regions but beyond this point, we find that the extent 
of the H II region increases as the density of the environ¬ 
ment decreases. After 10 Myr, we find that the ionization 
front has exited the box and all three regions display similar 
median neutral fractions. 

In the left panel of Figure we plot the positions 
and masses of the haloes found in a slice of our simula¬ 


tion the same thickness as the slice over which we project 
the neutral fraction. Here we use a simulation of our over- 
dense region with initial volume-weighted neutral fraction 
(/h i)init = 0.02 at tq = 1 Myr, as the ionization front is 
still contained within our box at this time and for this ini¬ 
tial neutral fraction we see a large scatter in the distance the 
front has travelled. We find that sometimes the progression 
of the ionization front appears to be impeded by the presence 
of a halo. This is further investigated in the middle panel of 
Figure]^ where we plot the H i column density found at the 
edge of the ionization front against the distance that the 
front has travelled, coloured by the average neutral fraction. 
We define the edge of the front to be where the smoothed 
neutral fraction becomes larger than log /hi = —4. We find 
that the absorbers slowing the progress of the ionization 
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front tend to be set back from the edge of the front and that 
we need to integrate over a wide distance of ±500 km s“^ 
(or ±0.6 pMpc, ignoring the contribution of peculiar veloci¬ 
ties) around the edge of the ionization front to capture their 
contribution to the column density. These absorbers have 
H I column densities in the range log A^h i(cm“^) = 16 — 21. 
The self-shielding prescription results in a smooth increase 
of the neutral fraction of the absorbers. 

We have also looked at the total hydrogen column den¬ 
sities of these absorbers and find that the absorbers with 
the highest column densities tend to be in the central 2 
pMpc of the region, surrounding the host halo. Outside 
this region, we find hydrogen column densities in the range 
logiVH(cm“^) = 20.5 — 21. At a given distance from the 
QSO the absorbers with the highest neutral fraction have 
the highest total column densities, as expected. For com¬ 
parison, we plot the column density we would expect in our 
interval at the mean density for z = 7, which brackets the 
lower end of the scatter in column density that we measure. 

To investigate the mass of the dark-matter haloes as¬ 
sociated with the edge of the ionization front, we found the 
halo closest to the point along our sightline where we define 
the ionization front to end. We hnd that the closest dark 
matter haloes are preferentially of low-mass (< 10 ^‘^Mq, 
left panel of Figure [^. The distribution thereby roughly 
tracks the distribution of haloes as a function of mass found 
in our simulation with a moderate enhancement at the low- 
mass end. There are also a significant number of sight lines 
whose closest halo is the QSO host halo. These correspond 
to directions in which the ionization front is still trapped 
within the host halo after tg = 1 Myr. 

In the middle panel of Figure we examine how the 
distance from the edge of the front to the closest halo, d, 
changes as a function of halo mass. We normalise this dis¬ 
tance by the virial radius i? 2 oo, which we define as the point 
where the density is equal to 200 times the mean enclosed 
density. We also choose random points within our simula¬ 
tion and found the closest dark matter halo. In this case, 
the median distance to the closest halo was slightly larger 
than this distance from the haloes to the ionization front, 
but we hnd a much larger scatter in the measured distances 
to the nearest halo. There is therefore only a rather loose 
but still signihcant association of the location of the ion¬ 
ization front and dark matter haloes, and this will depend 
both on where we dehne the point where the ionization front 
end and how we choose the halo that is associated with the 
edge of the ionization front. We have also checked how the 
H I and total column density of the absorber scales with the 
distance to the nearest halo d (right panel of Figure]^, but 
find no strong trend for either quantity. 


3.2 Synthetic Spectra 

We used our simulation to produce mock spectra, from 
which we measured the QSO near-zone sizes and compared 
to the observed near-zones. The mock Lyman-a spectra 
were constructed with rest-frame wavelength Aa = 1215.67 
A, damping constant Ka = 6.265 x 10® s“^ and oscillator 
strength fa = 0.4164. The peculiar velocity and tempera¬ 
ture of the gas were taken into account. To make our spectra 
more realistic, we added noise and convolved them with a 
Gaussian instrument profile with a FWHM of 240 km s“^. 


Figure shows examples of spectra constructed along 
three lines of sight in our over dense region at z = 7, vary¬ 
ing both the initial average volume-weighted neutral frac¬ 
tion and quasar lifetime. We find that, as expected from 
Figure the Lya transmission at fixed lifetime extends for 
longer distances in a more ionized IGM. The spectra con¬ 
structed from simulations with (/h i)init = 0.02 and 0.1 show 
broadly similar transmission inside the H il region, while 
the spectra from the (/h i)init = 1-0 display strong damp¬ 
ing wings. When the initial neutral fraction is kept con¬ 
stant at (/h l)init — 0.1 and the QSO lifetime varied (right 
panel of Figure]^, we see that along some sight lines (mid¬ 
dle and bottom panels) the front initially grows rapidly but 
the growth saturates already after 1 Myr. When there is a 
dense absorber along the sightline (top panel), saturation 
of the radius of the front is not yet reached even after 10 
Myr. In this case, the absorber has a total hydrogen column 
density log A^H(cm“^) = 21.5 (integrating over the same ve¬ 
locity interval as in Section 1^. The density begins to rise 
rapidly after a distance 3.4 pMpc from the source (which 
corresponds to the edge of the ionization front) before peak¬ 
ing at 3.7 pMpc. 

We determined the extent of the near-zone along each 
line of sight, which, following Fan et al. (2006), we defined 
as the distance from the source to the point where the trans¬ 
mitted flux drops below 10 per cent for the first time after 
smoothing our spectra with a box car window of width 20 
A in the observed frame. This smoothing reduces the effect 
of noise and intervening absorption lines. 


3.3 Near-Zones vs. Proximity Zones: Analytic 
Models 

It is important to note that the size of the H ll region is not 
necessarily the size of the near-zone that will be observed. 


Bolton & Haehnelt (2007a) summarise the analytical solu¬ 


tion for the radius of an ionization front moving into the 
IGM (i?ion). For an source emitting N ionizing photons per 
unit time, the rate of expansion of the H ll region is given 
by 


di?ion _ ^ - fTT-RionaH ii^h 
di lUH 


( 1 ) 


where /h i is the neutral hydrogen fraction, nu is the total 
hydrogen number density and an n is the recombination co¬ 
efficient for ionized hydrogen. Note that we do not use the 


case A recombination rate as in Bolton & Haehnelt (2007a), 


but the case B recombination rate as we are using the “on- 
the-spot” approximation for recombination. This difference 
is discussed further in Section [3^ Solving for i?ion, and as¬ 
suming that the lifetime of the source tq is much smaller 
than the timescale for recombination Gee = (^H<aHii)~^, 
leads to the solution 


T^ion — 


1.48 


N 


(A/hi)i/ 3 V1.3 X lO^^s-i 
X I —— 1 pMpe. 


1/3 


Vio«yr/ 


( 2 ) 


Due to residual neutral hydrogen inside the H ll region, 
the absorption can become saturated closer to the source. 
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Figure 6. Left: Simulated spectra through the overdeuse regiou for iuitial ueutral fractiou (/h i}init = (1.0,0.1,0.02) at tg = 1 Myr. 
Right: Simulated spectra for iuitial ueutral fractiou {/h i}init = 0.1 at tg = (1, 3, 10) Myr. Each row displays a spectrum aloug a differeut 
sightliue. The horizoutal dotted hue marks the poiut where the trausmitted flux is equal to 10 per ceut aud the vertical dotted hue marks 
the positiou of rest-frame Lymau-a. The spectra have beeu rebiuued outo pixels width 0.4 A iu the rest frame. 


By choosing the /h i that corresponds to the limiting opti¬ 
cal depth and assuming that the hydrogen gas beyond the 
ionization front is in ionization equilibrium, the largest ob¬ 
servable near-zone size, can be defined by 


T^max 


1.86 

^lim 


X 


f Tlim 


N 


1.3 X lO^^s-i 


1/2 


\ C 


1/2 


OL ^ [ci + 3] 


2 X lO^K ) 

- 1/2 


■(T 


-9/4 


pMpc, 

(3) 


where Aum is the overdensity corresponding to the limiting 
neutral fraction of hydrogen, T is the temperature, rum is 
the optical depth detection limit and a is the spectral index 
of the QSO. Note that this scales as rather than 
as before and that there is no dependence on /h i or tg. We 
confirmed this scaling using our radiative transfer simula¬ 
tions for luminosities varying over a range of 1 dex below 
our fiducial luminosity. 


3.4 The Sizes of Simulated Near-Zones 

We investigate the dependence of the near-zone size on the 
mass of the host halo by following the evolution of the me¬ 
dian near-zone size with time in our three different environ¬ 
ments at z = 6 and 7 (Figure [^. Note that from this point 
on, due to the computational costs running the radiative 
transfer simulations to 10 Myr, we only perform the radia¬ 
tive transfer on one octant of our simulation box. We find 
that once the ionization front has moved beyond the point 
where the clustering around host halo boosts the density of 
the IGM, there is only a very weak trend of decreasing near¬ 
zone size with increasing mass of the host halo. The other 
point of interest is the manner in which the near-zone sizes 
grow. There is a period of initial rapid growth, followed by 
a levelling-off in which there is little evolution in the median 
near-zone size of the system. The period of rapid growth is 


the time for which Rion < R^^^- After this point, we find 
that the median near-zone size and scatter in the near-zone 
sizes are independent of the initial neutral fraction. Due to 
the size of our simulation, at z = 6 approximately 10 per 
cent of our sightlines have near-zones that extend beyond 
the size of our box. In this case, we take the size of the box 
to be a lower limit to the near-zone size. In Appendix [B] we 
investigate the effect of including thermal feedback (from 
AGN and supernovae) in our hydrodynamical simulations 
and find that although this does not change our results sig¬ 
nificantly for tg > 1 Myr, we do find larger near-zone sizes 
for QSO lifetimes shorter than this due to a bubble of hot 
gas centred on the AGN. We find also that there is a natural 
decrease in the distance at which the near-zone sizes become 
saturated with increasing redshift. As shown by equation 
such a decrease is expected due to the increase of the mean 
baryonic density with redshift, as Aum and the IGM temper¬ 
ature are probably only varying slowly. We see also a small 
decrease in the scatter with increasing redshift. 

In Figure we plot the near-zone sizes we measure 
after tg = 1,3 and 10 Myr as a function of the mass of 
the host halo at z = 7. Results are shown for a range of 
initial average neutral fractions. Again, we find only a weak 
correlation between the near-zone size and the density of the 
three regions. This can be easily understood by looking at 
the left panel of Figure which shows the mean overdensity 
contained in a sphere of a given radius centred around the 
host halo. When the size of the observed near-zones are small 
enough, the matter overdensity around the host halo has 
not yet returned to the cosmic mean and the effect of the 
immediate environment of the host halo is imprinted in the 
distribution of near-zone sizes. However, due to the large 
scatter in near-zone sizes found in all three cases, it will be 
diflflcult to use near-zone sizes to constrain the mass of the 
host halo of the QSO. 

For our overdense halo, we also investigate how our 
choice of self-shielding model effects the sizes of the near- 
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QSO Lifetime (Myr) QSO Lifetime (Myr) QSO Lifetime (Myr) 


Figure 7. Evolution of QSO near-zone size with time for the three different density fields. Left: average, middle: intermediate, right: 
overdense. Top: 2 : = 6, bottom: 2 : = 7. In each panel we plot a range of different initial neutral fractions. The solid line in each case is 
the median near-zone size we measure for the system. The shaded region in the top panels shows observed near-zone sizes, rescaled to 
match the luminosity of our source, in the range 5.9 < 2 : < 6.1 ( [Carilli et al.|2010| > and in the bottom panels it shows the estimate for 
the near-zone size of the 2 : = 7.085 QSO ULAS J1120+0641. 



log Mhaio (h" Mo) log (h-’ Mo) log (h" Mo) 


Figure 8. Median near-zone size as a function of mass of the host halo at 2 ; = 7. The shift in the points around the masses of the three 
haloes is for display purposes only. The dashed grey vertical lines mark the mass of the host halo. The error bars represent the 15^^/85^^ 
percentiles in the near-zone sizes. No strong trend of near-zone size with the mass of the host halo is visible. The shaded region shows 
the estimate for the near-zone size of the 2 ; = 7.085 QSO ULAS J1120+0641. 


zones we measure. For the same photoionization rate log Th i at slowing the progression of the ionization front. This re- 


(s-i) = -14.3, 

we compare our fiducial model, taken from 

suits in a smaller H ll region and near-zone sizes. However, 

Rahmati et al. 

(2013), with a simple threshold self-shielding 

the distance at which the near-zone sizes saturate is similar 

model where all gas with overdensity above some thresh- 

for both self-shielding models and for the case of no self¬ 

old Ass (|Schaye 2001) is assumed to be neutral (see also 

shielding at all. Note also that the largest effect of chang¬ 

fMiralda-Escude, Haehnelt & Rees||2000| |Furlanetto & Oh| 

ing the self-shielding model is the amount of neutral gas 

120051 Bolton et aL|2011 ). We find that in the simple thresh- 

available for the same photoionization rate. The difference 

old model, the Lyman limit systems become more effective 

between the threshold self-shielding model and the prescrip- 
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tion of Rahmati et al. (2013) increases the average volume 


weighted neutral hydrogen fraction by an order of magnitude 
and neglecting self-shielding entirely decreases the average 
volume weighted neutral hydrogen fraction by an order of 
magnitude. 

We also check the effect of switching from case B to case 
A recombination rates. We find that the near-zone sizes sat¬ 
urate at a distance of 3 pMpc from the host halo, closer 


(2013) model. We also do not find a large change in the 


to the value obtained in Bolton et al. (2011) than what we 


find in our models using case B recombination. However, 
we still do not include the diffuse radiation that would re¬ 
sult from a full treatment of recombinations. The change 
in recombination rate does not make a large difference for 
the smaller near-zones we measure and only becomes impor¬ 
tant along sightlines where the ionization front has travelled 
a large distance when the number of recombinations along 
the sightline becomes significant. This therefore also has the 
result of reducing the scatter slightly. 

3.5 Red Damping Wings in the Simulated Spectra 

The second observable that has been used to quantify the 
impact of a still largely neutral IGM on the spectra of 
QSOs is the transmission at the rest-frame wavelength of 
Lya. In the presence of high densities of neutral hydrogen, 
even flux redward of Lya can be attenuated producing a red 
damping wing. The strength of this damping wing can then 
be used to place constraints on the neutral fraction of the 


IGM ( Miralda-Escude|1998 Mesinger, Haiman Gen|2004 


Wyithe Loeb||2 004| [Me singer Furlanetto||20Q8| Bolton 


et al.||201H [Mortlock et al.||201l] [Schroeder, Mesinger fc 

Haiman|2013|). 

In Figure we plot the transmission at the Lya rest- 
frame wavelength as a function of time for our three regions. 
We measure this transmission after smoothing our spectra to 
the same resolution at which we measure the near-zone sizes. 
We find that, as expected, the regions with higher average 
initial neutral fractions display stronger red damping wings. 
These damping wings are stronger for shorter QSO lifetimes, 
when the ionization front is still in the higher-density region 
surrounding the host halo. 

Perhaps more interestingly, we also find that the trans¬ 
mission at Lya is consistently lower in the overdense regions, 
for all the initial average neutral fractions that we consider 
due to the increased number of dense, self-shielded systems. 
This is displayed more clearly in Figure |10| It is perhaps 
not surprising that we find that the mass of the host halo 
plays an important role in the strength of the red damping 
wing and not for the size of the near-zones. In the case of 
the near-zones, we are probing parts of the IGM with aver¬ 
age density close to the mean. The strength of the damping 
wing is however more sensitive to the environment around 
the source. This clustering of dense absorbers is largest in 
the most massive halo we look at (left panel of Figure |^, 
resulting in a large scatter in the transmission at at Lya not 
found in the lower mass haloes we simulate. 

In this case, the choice of self-shielding model does not 
seem to make a large difference. For short QSO lifetimes, 
we find that the threshold self-shielding model predicts a 
stronger damping wing. For longer, more realistic lifetimes 

10 Myr), the results from the threshold self-shielding 
model are similar to those predicted by the [Rahmati et ah] 


transmission when we use the case A recombination rate. 


4 COMPARISON WITH OBSERVATIONS 
4.1 z ~ 6 Near-Zones 

In the left panel of Figure m we compare the measured 
near-zone sizes from simulation of our overdense region at 
z = 6 to observed near-zone sizes in the range 5.9 < z < 6.1 


taken from Garilli et al. (2010). Here, we used a background 
photoionization rate logFHi(s“^) = —12.8, based on mea¬ 
surements of the UV background at z ^ 6 by jGalverley et aQ 
(2011) and Bolton & Haehnelt (2007a) and logFHi(s“ ) = 
— 12.5 to mimic a late reionization scenario (as in [Ghoud- 


hury et ^|2014| . We consider two different QSO luminosi¬ 


ties, log (s ^) = 57.1 and 57.3, which fall in the range of 
values estimated by Fan et al. (2006). 


For the same quasar luminosity and background pho- 
toionization rate, there is an increase in the near-zone size 
with decreasing redshift. With our fiducial quasar luminos¬ 
ity and a UV background logFni = —12.8, we match the 
lower end of the observed z = 6 near-zone sizes remarkably 


well. This is somewhat in contrast with other work (Bolton 
& Haehnelt 2007a Wyithe, Bolton V Haehnelt||20(38 ) and 


the different result is most likely due to the difference in the 
recombination rates used here (case B vs. case A). We have 
also run a simulation where we increased the luminosity of 
our ionizing source by a factor of 1.5 and find that this in¬ 
creases the near-zone sizes we find with many now larger 
than the near-zone sizes our box can contain. In the cases 
where the near-zone size occurred at a distance larger than 
the radius of our box, we took the radius of the box as a 
lower limit. 


4.2 The z = 7.085 QSO ULAS J1120+0641 

We now want to use our simulations to revisit the constraints 
on the ionization state of the IGM around the z = 7.085 
QSO ULAS J1120+0641. Modelling of this kind has pre¬ 
viously been undertaken in Bolton et al. (2011). However, 


the mass of the host halo in that work was a factor of 100 
times less massive than the one considered here (although 
based on section [3^ we expect this to make remarkably little 
difference). 

Our analytical estimate described by equation © pre¬ 
dicts that = 1.86/Alim proper Mpc (assuming our 

fiducial values for a and N). Assuming that the absorption 
at z = 7.085 is due to overdensities of A ~ 0.5 similar to that 
required to explain the observed near-zone sizes at z = 6, 
then we would expect ~ 3.72 (albeit with a large scat¬ 


ter). Mortlock et al. (2011) measure the near-zone of ULAS 
J1120+0641 to be 1.9 ± 0.1 proper Mpc. This suggests that 
the near-zone of ULAS J1120+0641 is unlikely to be due to 
the proximity effect and may therefore be a useful probe of 
a substantially neutral IGM at z = 7.085. 

In the right panel of Figure we compare the mea¬ 
sured near-zone size of ULAS J1120+0641 (shaded purple 
region) to the simulations of the overdense region containing 
a 2.5 X 10^^ Mq halo. We find that only for simulations with 
a large average initial neutral fraction somewhere between 
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Figure 9. Evolution of transmission at Lya with time for the three different density fields. Left: average, middle: intermediate, right: 
overdense. Top: 2 : = 6, bottom: 2 : = 7. In each panel we plot a range of different initial neutral fractions. The solid line in each case is the 
median transmission at 1216 A we measure. More sightlines display a red damping wing in the overdense region than in the less dense 
regions. The shaded region shows the estimate for the transmission at Lyo; of the 2 : = 7.085 QSO ULAS J1120+0641. 



9 10 11 12 13 9 10 11 12 13 9 10 11 12 13 
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Figure 10. Median transmission at Lyo: as a function of mass of the host halo at 2 : = 7. The shift in the points around the masses of 
the three haloes is for display purposes only. The dashed grey vertical lines mark the mass of the host halo. The error bars represent the 
15^^/85^^ percentile scatter in the near-zone sizes. The shaded region shows the estimate for the transmission at Lyo; of the 2 : = 7.085 
QSO ULAS J1120+0641. 


(/h i) init = 0.1 — 1.0 can we reproduce the small observed 
near-zone size. 

For QSO lifetimes long enough to allow the central mas¬ 
sive black hole to grow while the QSO is optically bright, 
we also struggle to match the strength of the damping wing 
observed in ULAS J1120+0641 except in the case of our 
most massive host halo or for an initially neutral IGM. For 
the simulations centred on our most massive halo, damping 
wings almost as strong as that of ULAS J1120+0641 are 


seen in the 15^^/ 85^^ percentiles of all of our models, even 
the ones with low volume weighted initial neutral hydrogen 
fraction (e.g., (/Hi)init = 0.02). The observed red damping 
wing therefore provides little constraint on the neutral frac¬ 
tion of the surrounding IGM if ULAS J1120+0641 is indeed 
hosted by a massive halo. However, note here that we have 
neglected the ionizing radiation from Lyman-a emitters that 
may reside in these systems and weaken the Lya damping 
wing ( Bolton fc Haehn^|2013 ). 
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Figure 11. Left: QSO near-zone sizes (corrected for differences in luminosity) at 2 : = 6 and z = 7 for background photoionization rates 
(logPH I (s“^) = -12.8, -12.5) and two different QSO luminosities (log A^-^(s“^) = 57.1, 57.3) in our overdense region after 10 Myr. The 
shaded region represents the observed near-zone sizes at 2 ; = 6 taken from [Carilli et aL] ( |2010| . Middle: Observed near-zone sizes as a 
function of redshift, scaled to the same luminosity. Plotted in blue is the analytic evolution of near-zone sizes with redshift given by 
Equation with Aii^ at the median and the measured 15^^/ 85^^ percentiles of the saturated near-zone sizes at 2 ; = 7. The shaded 
regions are the 15^^/85^^ percentiles of the simulated near-zone sizes at 2 ; = 6, 7. Right: Distribution of the near-zone sizes at 2 ; = 6 and 
7 in the simulations with logEn i (s“^) = -12.8 after 10 Myr for our default luminosity. 


Note also that modelling of this kind is further com¬ 
plicated by the fact that, towards the end of reionization, 
there are expected to be large variances in the neutral frac¬ 
tion along different sightlines (e.g., Becker et al.||2015 ) and 
that the QSOs may have been bright at earlier redshifts and 
may have already ionized its surroundings ( Feng et al.|2013 ). 


4.3 Evolution with Redshift 


In the middle panel of Figure mi we plot the published 
near-zone sizes of a selection of quasars with z > 5.8 
([Carilli et al. 2010[ [Willott et al.||2010[ [Mortlock et al. 


2011 


|Reed et al.||2015| [Venemans et al.||2015| |Wu et al. 


2015). Near-zone sizes have been corrected for luminosity 
as flNz,corr = 10° ‘‘^^^+^i450,AB)/2. Note that near-zone 
sizes in the literature have been often corrected assuming 
that the near-zone size scales with the ionizing luminosity 
as ]Vl/3 

, as expected when the near-zone size traces the edge 
of the H II region, not which is the (more likely correct) 
scaling expected for saturated near-zone sizes. The recently 
reported rather large near-zone size measured by |Wu et al.| 
(2015) decreases substantially when rescaled for the lumi¬ 


nosity of this very bright QSO. We plot this particular near¬ 
zone size as an upper-limit, as Wu et al. (2015) have used 


a definition of near-zone size for their measurement that is 
different from the normally used definition fromlFan et al.l 
( 200 ^ . 

For comparison, we also plot the values at the 15^^/85^^ 
percentiles of near-zone sizes we measure in our simulations 
at z = 6 and 7, also rescaled for luminosity (with an as¬ 
sumed Mi 45 o,ab to match ULAS J1120+0641). Here we use 
a highly ionized IGM (with background photo-ionization 
rate logFn i (s“^) = -12.8, to try and make a fair compari¬ 
son to the observed z — e near-zones). For this Fh i, we find 
that the near-zone sizes saturate quickly. Here we plot the 
scatter of a simulation at tq = 10 Myr, but in this regime 
we no longer expect any dependence on tq (Figure]^. We 
also plot the analytic solution for the evolution of saturated 


near-zone sizes with redshift (Equation]^, where we have 
selected Aum to fit the median near-zone size and the scat¬ 
ter we measure at z = 7 [Aum = (0.26, 0.38, 0.48) for the 
15^^ percentile, median and the 85^^ percentile]. 


The analytic solution tuned to fit our results at 2 ; = 7 
also agrees well with our results at z = 6. Note that the 
evolution we predict is shallower than the linear fits to the 


near-zone sizes from Carilli et al. (2010) and Venemans et al. 


(2015) perhaps suggesting that the IGM indeed becomes 
substantially more neutral at z > 6.3 than predicted by the 
evolution of the mean density alone. We also note that the 
scatter at z = 6 will in reality be somewhat larger (indicated 
by the grey arrow in the plot), as we have used the size of 
our box as a lower limit in cases where we do not find a 
near-zone size along a sight line. At z > 6.5 there are several 
near-zones that fall outside of our predicted scatter. These 
near-zones may represent unsaturated cases, where the size 
is still following the radius of the H ll region. For realistic 
QSO lifetimes tq ^ 1 Myr, Figure [^suggests that these QSO 
may reside in regions of the IGM with (/Hi)init ^ 0.1 and 
lower photo-ionization rate than we assumed in Figure 


The observations at z ^ 6 contain near-zones with sizes 
>10 pMpc which our simulations cannot reproduce due to 
the limited box size of our simulation. The distribution of 
near-zone sizes at z = 6 shown in the right panel of Figure 
|TT] has a large peak at 8 pMpc. This is the result of the 
lower limits we are assuming for sightlines where no near¬ 
zone is found. For a larger simulated region some of our 
sightlines should produce near-zones that would match the 
observations. However, these would be the outliers in our 
distribution and it is also possible that defining the point 
where the near-zone ends at z ~ 6 becomes ambiguous due 
to patches of transmitted flux in the Lya forest (Bolton & 
Haehnelt]|2007a ). 
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4.4 Caveats and Limitations 


There are obviously still a significant number of uncertain¬ 
ties and limitations in our modelling that should be taken 
into account when considering our results. First, the reso¬ 
lution of our hydrodynamical simulations is still somewhat 
below what would ideally required to resolve the Lyo forest 


(Bolton & Becker 2009) at 2 : > 5. Our maximum resolu¬ 


tion was dictated by the large simulation region necessary 
to follow the evolution of the ionization front over a realis¬ 
tic quasar lifetime. For the simulation box to still contain 
observed near-zone sizes at z > 6.5 , the chosen resolution 
was the best that was still computationally feasible to run. 
As we have discussed, the largest observed near zone sizes 
at z ~ 6 still exceed the size of the region we have simulated 
here by up to a factor 1.5. We have investigated the effects 
of changing the resolution in Appendix [A] 


As stated in Section |2.2[ our radiative transfer post¬ 
processing only evolves the temperature and ionization state 
of the gas, but neglects the evolution of the density field as 
well as the coupling of the radiation to the density field. 
However, as discussed earlier the timescales that we are 
considering here are short compared to the timescales over 
which this matters. The speed of the ionization front in our 
case is close to the speed of light, and therefore many or¬ 
ders of magnitude larger than the sound speed in the IGM. 
Another possible effect we neglected is the thermal pressure 
due to the heating provided by the passage of the ionization 
front. But, once again, adiabatic expansion of the gas is 
negligible over the timescales that we are considering. Also, 
as we are using the “on-the-spot” treatment of recombina¬ 
tion radiation, we assume that ionizing photons produced 
by hydrogen and helium recombinations are re-absorbed in¬ 
stantaneously in the same computational cell. The inclusion 
of the diffuse radiation from recombinations would likely re¬ 
duce the ability of dense absorbers to effectively shadow the 
ionization front (as seen in Figure]^ and somewhat reduce 
the scatter in the near-zone sizes we model [see Cantalupo 
& Porciani ( 2Q11| for a detailed discussion of recombination 
radiation effects on quasar H II regions]. 

We furthermore only model the initial ionization state 
of the IGM using a homogeneous UV background. This 
certainly oversimplifies the previous evolution of the IGM 
before the QSO turns on for a single episode of constant 
luminosity as reionization is a highly inhomogeneous pro¬ 
cess, with individual galactic sources pre-ionizing inhomoge- 
nously the surrounding IGM to different levels and the QSO 
likely to have a complex light curve. An enhanced contribu¬ 
tion to the photoionization rate from galaxies close to the 
QSO has been previously suggested to increase near-zone 
sizes at z = 6 (Bolton & Haehnelt 2007a| Wyithe, Bolton 


fc Haehnelt||2QQ8 Maselli, Ferrara Gallerani[|2009 ). How¬ 

ever, this simple model for the UV background should not 
change our predictions for the evolution of saturated near¬ 


zone sizes with redshift (middle panel of Figure 11), as the 
distance at which the near-zone sizes become saturated at 
is independent of the ionization state of the gas. It will also 
not change our results at z = 7, as we find it difficult to 
match the observed near-zone of ULAS J1120+0641 even 
for an IGM that is initially completely neutral. 


5 SUMMARY AND CONCLUSIONS 

We have presented here results from our detailed modelling 
of quasar near-zones at z = 6 and 7. We have performed 
zoom simulations of three dark matter haloes from the Mil¬ 
lennium simulation, which differ in mass by more than a 
factor of 100, and post-processed these with the 3D radia¬ 
tive transfer code RADAMESH, exploring a range of initial 
neutral fractions representing different reionization histories. 
We placed a source of the same luminosity as expected for 
the = 7.085 QSO ULAS J1120+0641 in each halo. 

The spatial extent of the ionization fronts in our sim¬ 
ulations shows a large scatter with direction for fixed QSO 
lifetime. The location of the edge of ionization front appears 
to have a loose association with low-mass dark matter haloes 
(> IQpM q) hosting self-shielded absorbers. We further find 
that the size of the near-zones of bright high-redshift QSOs 
depends only very weakly on the mass of the QSO host halo. 
For realistic QSO lifetimes (> 1 Myr), the ionization front 
has moved out into the general IGM where the average den¬ 
sities have returned to close to the mean. 

Our simulations reproduce the observed range of z ~ 

6 — 7 near-zone sizes very well once we take into account 
the limitations due the size of our simulated region. Even 
when the extent of the near-zone sizes in our simulation has 
become saturated due to residual neutral hydrogen in the 
H II region, we still hnd a large scatter in the near-zone sizes 
along different sight lines around the host halo. The scatter 
of the observed flux due to variations in the density held 
can explain a lot of the scatter in the observed near-zone 
sizes without the need to invoke different volume averaged 
neutral hydrogen fractions along different lines of sight or 
anisotropic emission. 

We have presented here the hrst simulations of QSO 
near-zones around a host halo with a mass as high as ex¬ 
pected for ULAS J1120+0641. We hnd that, for a QSO 
lifetime of 10 Myr, the observed near-zone size of ULAS 
J1120+0641 falls just within the lower end of our 15^^/85^^ 
percentile scatter for a IGM that is initially completely neu¬ 
tral at z = 7. This may be indicative of the presence of a 
high H I column density absorber lying close to the quasar or 
may suggest that the duration of the optically bright phase 
in ULAS J1120+0641 is < 10 Myr and thus too short for 
its supermassive black hole to grow signihcantly. Alterna¬ 
tively, this may cast doubt on the robustness of the mea¬ 
sured size of this particular near-zone. We further caution 
that for massive QSO host haloes, the transmission at the 
Lya rest-frame wavelength can be low even if the IGM is 
generally highly ionized. 

Gurrent and upcoming surveys will provide new obser¬ 
vations of QSOs and their near-zone sizes which, in conjunc¬ 
tion with alternative observations, should allow us to further 
improve constraints on the neutral fraction of the IGM at 
the end of reionization, as well as on the relative contribu¬ 
tion of obscured and unobscured growth to the build-up of 
the supermassive black holes powering them. 
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APPENDIX A: RESOLUTION TESTS 


To study the effect of the numerical resolution of our simu¬ 
lation on the near-zone sizes we measure, we ran two addi¬ 
tional hydrodynamic simulations of our overdense regioiQ 
differing in mass resolution from our fiducial model by a fac¬ 
tor of eight, summarised in Table [AT] For the high resolu¬ 


tion run, only a sphere with radius 9 comoving Mpc was 
simulated due to computational restraints. The high (low) 
resolution hydrodynamical simulations were mapped onto 
correspondingly finer (coarser) grids with respect to our de¬ 
fault model. Examples of the resulting ionization fronts in a 
simulation with (/h i)init = 0.02 at tq = 0.2 Myr are shown 
in Figure |A1[ We find that the scatter in the distance the 
ionization front has travelled among the different sightlines 
increases as the resolution is increased. 

Figure shows the results of three tests we performed 
to understand how changing the resolution altered our mea¬ 
surements of the near-zones and the transmission at the Lya 
rest frame wavelength. In the left panel, we show the median 
near-zone size and the values at the 15^^/85^^ percentiles 
as a function of the gas particle mass. The middle panel 
shows the frequency of different near-zone sizes among the 
sight lines. We use a snapshot at tq = 2 and 10 Myr in 
a simulation with (/Hi)init = 1-0. For the simulation with 
tq = 10 Myr, we only plot the near-zones measured in our 
default and low resolution simulations, as for this QSO life¬ 
time many of the near-zone sizes are larger than the radius of 
our highest resolution simulation. We find that, as the reso¬ 
lution is degraded, our simulations predict a similar median 
and scatter in near-zone sizes. 

We further find that the transmission at Lya is similar 
in our default and low resolution simulations at both QSO 


lifetimes considered here (right panel of A2). The transmis¬ 
sion at the Lya rest-frame wavelength is higher in our high¬ 
est resolution simulation, most likely due to our star forma¬ 
tion prescription. As we increase the resolution of our simu¬ 
lations, we form stars more efficiently and there is thus less 
gas remaining. This will have the largest effect in the dens¬ 
est regions responsible for absorption. These are the regions 
that also should be affected by stellar and AGN feedback 
not included in our default simulations. We note, however, 
that the median transmission and the scatter in our default 
simulation is reasonably well matched by a simulation with a 
more realistic star formation prescription that also includes 
stellar and AGN feedback (Section [b|. 


APPENDIX B: EFFECT OF FEEDBACK 

As we expect feedback from stars and AGN to play a role 
in shaping the environment around galaxies, we ran another 


^ In each of these appendices, we perform tests using our over- 
dense region at 2 ; = 7 which contains the 10^^-^ Mq halo, unless 
stated otherwise. 


simulation of our overdense region which included the ef¬ 
fect of feedback. We have used here the same feedback pre¬ 
scription as in Gosta, Sijacki & Haehnelt (2015), including 
supernovae- and AGN-driven outflows, which reproduced 
well observations of spatially extended cold gas around a 
2 : = 6.4 QSO. For the supernovae-driven winds, we used a 
mass loading factor rj — 1 and a velocity of 1183 km s“^. 
However, note that our simulation fails to reproduce the fast- 
moving cold gas found in Gosta, Sijacki & Haehnelt (2015), 
perhaps due to differences in the codes used (the moving- 
mesh code AREPO versus the SPH code p-GADGEt3), a dif¬ 
ferent choice of UV background or the effect of including 
metal-line cooling. 

We find that the shape of the ionization front is re¬ 
markably similar to our fiducial model (compare the middle 
and right panels of Figure Bl). There is still a large scat¬ 
ter in the distance that the front has travelled, depending 
on the sight line. In the middle panel of Figure [B^ we plot 
the growth of the median near-zone size with time. We find 
that, for tq < 0.5 Myr, the model including feedback shows 
near-zone sizes with median size ~ 0.5 pMpc larger than our 
fiducial model. This is due to a combination of the galactic 
winds removing dense gas from the host halo and the ther¬ 
mal feedback resulting in more collisionally ionized gas in 
the host halo (also visible in some of the larger haloes in the 
left panel of Figure [bT] ). The right panel of Figure B2 shows 
that the temperature of the gas surrounding the host halo is 
an order of magnitude higher in the case with the AGN and 
SN feedback compared with the simulation run with Quick 
Lyman-a star formation. 

However, does not have a large effect on the Lya trans¬ 
mission and for longer QSO lifetimes we find that the median 
near-zone size and scatter found in our fiducial model is very 
similar to that in the simulations with feedback. The star 
formation prescription we have used in our fiducial model 
prevents the build-up of high-density gas in the host halo. 
As well as being computationally efficient, this has the ef¬ 
fect of mimicking more physically motivated feedback mod¬ 
els which will eject gas from the host halo. We also find 
evidence for strong damping wings in the simulation includ¬ 
ing feedback (and note that the H i photoionization rate 
used here results in a volume weighted neutral fraction of 
0.02 in our fiducial model). 


APPENDIX C: LY/3 NEAR-ZONES 

We have also investigated the dependence of Ly/3 near-zones 
on the environment surrounding the host halo. [Bolton ^ 
Haehnelt ( 2007a|b ) have highlighted the potential of Lyj3 
absorption spectra as a probe of the z > 6 IGM, as they 
are expected to take longer to reach than Isya near¬ 

zones. The corresponding Lya forest absorption will occur 
at z' = [A^(l + z)/\a\ — 1. For z = 7, this corresp onds 
to absorption at z — 5.75. From Becker et al. (2015), the 


effective optical depth at this redshift will be greater than 
2, with some observations returning a lower limit of 7. To 
test the effect of changing the effective optical depth, we 
simply increase the Ly/3 optical depth by a constant To,{z') 
across the whole sightline. Note however taking Tc{z') to be 
constant is an unrealistic assumption as we expect Ta{z) to 
vary with changes in the density field, etc. In Figure [CT] we 
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X [pMpc] X [pMpc] X [pMpc] 


Figure Al. Map of the projected neutral hydrogen fraction over a 0.03 pMpc slice in our overdense region at tg = 0.2 Myr at three 
different resolutions in a simulation with {/h i}init = 0.02. Left: mgas = 2.2 x 10^ ^o- Middle: mgas = 2.8 x 10® h~^ Mq (default 

resolution). Right: mgas = 3.5 x 10® h~^ Mq. 



^gas(^o*^ ) Near-Zone Size (proper Mpc) T^216 

Figure A2. Left: median near-zone size with error bars representing the values at the 15^^/85^^ percentiles against the mass of the gas 
particle in the three different resolution simulations we test at tg = 2 Myr and for our default and low resolution simulations at 10 Myr 
with {/h i}init = l-O- Middle: histogram of the near-zone sizes measured in our overdense region with{/H i}init = 1-0 at tg = 2 and 10 
Myr. Right: histogram of the transmission at the Lyo; rest-frame wavelength measured in our overdense region with (/h i}init = 1-0 at 
tg = 2 Myr and 10 Myr. 


plot the Lya and Ly/3 spectra along a sightline through our 
simulation with ra{z') = (0,2,4). 


We find that, due to the smaller scattering cross-section 
of Ly^^, we see more transmitted flux in the Ly^^ spectrum 
than the Lya spectrum when Ta{z') = 0. This results in a 
correspondingly larger near-zone. When we tested the ratio 
of Ly/3 near-zone sizes to Lya near-zone sizes, we reproduce 
the results of Bolton & Haehnelt (2007b), where a smaller 
volume weighted neutral hydrogen fraction gives a larger 
Ly^/Lya near-zone ratio. 


However, when we include the additional contribution 
to the optical depth due to lower redshift Lya absorption, 
we find that the flux drops accordingly. We still find some in 
the Ly/3 region for ra{z') = 2, but for Ta{z') = 4 we see no 
significant transmission of flux along this sightline and Ly/3 
near-zones lose their effectiveness as a probe of the neutral 
fraction at z = 7. They may, however, still act a a valuable 
probe at lower redshift where the “contamination” of the 
Ly/3 region by Lya absorption is lower. 


APPENDIX D: HELIUM IONIZATION FRONTS 

We have also investigated how the helium ionization fronts 
progress and the effect that this as on the temperature of 
the IGM surrounding the QSO. For a simulation with a large 
photoionization rate T (here we present results for log Th i 
( s“^) = -13.7), we find that most of the helium is in He ll. 
Otherwise, for smaller T, we see a mixture of He I and He ll, 
with the He I found in the higher density regions. In the 
left panel of Figure [DT] we present a projection of the He ll 
fraction after 10 Myr. We find that the front has travelled 
a much smaller distance than say the front ionizing the H i 
(top panel of Figure]^ and the shape is reminiscent of the 
propagation of a H ll front into a neutral IGM (bottom panel 
of Figure]^. 

In the middle panel of Figure [DT] we show the median 
He II fraction across our sightlines at tg = 1, 3 and 10 Myr. 
Again, as with the case of the simulations with where the 
hydrogen in the IGM is initially neutral, we find that the 
front continues to grow steadily over the 10 Myr we simulate. 
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Figure Bl. Left: total gas mass contained in a sphere with radius specified on the horizontal axis, centred on the host halo for the 
two simulations. Middle: map of the projected neutral hydrogen fraction over a 0.03 pMpc slice for our overdense region for the Quick 
Lyman-o; star formation prescription at tg = 1 Myr with background photoionization rate loghHi (s“^) = -14.3. Right: map of the 
projected neutral hydrogen fraction for the simulation containing including AGN and SN feedback. 




Radius (proper Mpc) QSO Lifetime (Myr) QSO Lifetime (Myr) 

Figure B2. Left: mass-weighted average temperature in a sphere with radius specified on the horizontal axis, centred on the host halo for 
the two simulations. Middle: evolution of the near-zone sizes with time in our overdense region for the models with the Quick Lyman-o; 
star formation prescription and with AGN and SN feedback. The solid line represents the median value and the dashed lines are the 
values at the 15^^/85^^ percentiles. Right: evolution of the transmission at the Lyo rest frame wavelength with time. 
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Figure Dl. Left: projected map of the He ll fraction after Iq = 10 Myr in the simulation with an average initial neutral hydrogen 
fraction of {/h i}init = 0.002. Middle: the median neutral He ll fraction shown at tg = 1 Myr (black), Iq = 3 Myr (red) and tg = 10 
Myr (blue). Right: the median temperature along our sightlines for the simulation with (/h i}init = 0.002 shown at tg = 1 Myr (black), 
tQ = 3 Myr (red) and tg = 10 Myr (blue). 
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Resolution 

mgas(xl0®/i ^ Mq) 

Softening length (h ^ kpc) 

Size of largest cell (h ^ Mpc) 

High 

0.35 

0.625 

0.044 

Default 

2.80 

1.250 

0.089 

Low 

22.0 

2.50 

0.177 


Table Al. Comparison of the three simulations we consider in our resolution tests. 


Distance (proper Mpc) 



6.85 6.90 6.95 7.00 

Redshift 


Figure Cl. Simulated Lyo; and Ly/3 spectra along a sightline for 
initial neutral fraction {/h i}init = 0.02 at tg = 1 Myr. In each 
panel we consider a different value for the optical depth of the Lya 
forest at the redshift z' = [Xp{l z)/\a\ — 1 corresponding to 
the Ly/3 absorption. For Ta{z') = 0 (top panel), the contribution 
from the lower redshift Lya forest is neglected. For Toi.{z') = 2,4 
(middle and bottom panels), we experiment with different levels 
of contribution from the optical depth from the 2 ; = 5.75 Lya 
forest. Varying Tol {z') results in a dramatic change in the observed 
Ly/3 flux. 


At tq = 1 Myr, the front has travelled just over 1 pMpc 
from the host halo. This corresponds to the region of highly 
ionized hydrogen around the host halo (right panel of Figure 
1^ . We also find a low He ll fraction in the immediate vicinity 
of the host halo where the gas is very hot. 

We find that although there is a temperature boost as¬ 
sociated with the passage of the He ii front, it provides an 
increase of about 20 per cent only (right panel of Figure 
Dl). This will not have much of an effect on the recombi¬ 


nation rate and hence the near-zone sizes we measure. The 
radius of the thermal proximity zone is similar in the three 
regions around our three host haloes and we measure the 
same temperature at mean density around all three of our 
host haloes (where To is the mean volume-weighted gas tem¬ 
perature with a density within 5 per cent of the mean den¬ 
sity). At z = 6, we find that log To (K) = 4.30, within the la 
error bars of the measurement of the temperature around 
the z = 6 quasar SDSS J0818+1722 by [Bolton et al. (2010). 
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